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ABSTRACT 

The common approach to compute the cosmic-ray distribution in an starburst 
galaxy or region is equivalent to assume that at any point within that environment, 
there is an accelerator inputing cosmic rays at a reduced rate. This rate should be 
compatible with the overall volume-average injection, given by the total number of 
accelerators that were active during the starburst age. These assumptions seem rea- 
sonable, especially under the supposition of an homogeneous and isotropic distribution 
of accelerators. However, in this approach the temporal evolution of the superposed 
spectrum is not explicitly derived; rather, it is essentially assumed ab-initio. Here, we 
test the validity of this approach by following the temporal evolution and spatial dis- 
tribution of the superposed cosmic-ray spectrum and compare our results with those 
from theoretical models that treat the starburst region as a single source. In the 
calorimetric limit (with no cosmic-ray advection), homogeneity is reached (typically 
within 20%) across most of the starburst region. However, values of center-to-edge 
intensity ratios can amount to a factor of several. Differences between the common 
homogeneous assumption for the cosmic-ray distribution and our models are larger in 
the case of two-zone geometries, such as a central nucleus with a surrounding disc. We 
have also found that the decay of the cosmic-ray density following the duration of the 
starburst process is slow, and even approximately 1 Myr after the burst ends (for a 
gas density of 35 cm -3 ) it may still be within an order of magnitude of its peak value. 
Based on our simulations, it seems that the detection of a relatively hard spectrum up 
to the highest gamma-ray energies from nearby starburst galaxies favors a relatively 
small diffusion coefficient (i.e., long diffusion time) in the region where most of the 
emission originates. 

Key words: ISM: supernova remnants — cosmic rays — galaxies: starburst - 
gamma rays: galaxies — radiation mechanisms: non-thermal 



1 INTRODUCTION 

From radio to TeV gamma-rays, starbursts shine in the sky. 
Their large star formation rates power high injection rates 
of cosmic rays in the interstellar media of these galaxies, 
which emit nonthermal radiation. The large CR energy con- 
tents in starbursts are of interest for cosmic ray feedback 
(Socrates et al. 2008), ionization (Papadopoulos 2010), the 
gamma-ray background (e.g., Thompson, Quataert & Wax- 
man 2007; Makiya et al. 2011), and for determining whether 
they are neutrino sources (e.g., Loeb & Waxman 2006). In- 
deed, the large cosmic content has been supported by the 



recent detections of M82 and NGC 253 by GeV and TeV 
telescopes (Abdo et al. 2010, Acciari et al. 2009, Acero et 
al. 2009). In anticipation and as spinoff of these detections, 
there has recently been a spate of modelling of the cosmic 
ray spectrum in starburst galaxies. 

In considering the cosmic ray spectrum in star forming 
regions, models generally start by writing the diffusion-loss 
equation (see, e.g., Ginzburg & Syrovatskii 1964) 



-»v»w + 3g 



- A [b(E)N(E)] - Q(E) = 



dN(E) 
8^' 



(1) 
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where N(E) is the distribution of particles with energies 
in the range E and E + dE per unit volume, and b(E) = 
— (dE/dt) is the rate of loss of energy. In this equation, D 
is the scalar diffusion coefficient, Q(E) represents the source 
term appropriate to the production of particles with energy 
E, and t(E) is the residence time (beyond which particles 
are removed from the phase space). It is usually assumed 
that the starburst is in a steady state, dN(E)/dt — 0, 
and, under the assumption of a homogeneous distribution 
of sources, that the spatial dependence can be ignored, so 
that D v 2 N(E) = 0. These assumptions were adopted by 
essentially all detailed models of starburst regions published 
in the literature so far (e.g., see Paglione et al. 1996, Blom et 
al. 1999, Torres 2004, Domingo-Santamaria & Torres 2005, 
de Cea del Pozo et al. 2009, Lacki et al. 2010, and references 
therein). Persic et al. 2008, and Rephaeli et al. 2010 solve the 
diffusion-convection equation in 3D, but do not discuss the 
detailed spatial distribution of cosmic-rays, just the average 
spectrum over the starburst and the surrounding disk. This 
was also the case of some of the formerly quoted papers. 
On the other hand, the GALPROP models (see Strong & 
Moskalenko 1998) for the propagation of cosmic-rays in our 
Galaxy do not consider the formation of the cosmic-ray spec- 
trum from the individual contributors, in a time-dependent 
manner. 

In the previously mentioned models, and in regards of 
Eq. |T]), the proton injection is assumed to be a power-law 
with index p, for instance, 



tVinj (,-^p, 



= KE, 



p, kin 



(2) 



where K is a normalization constant and units are such that 



1= GeV" 1 cm" 



The proton injection is assumed to 



be uniform across the starburst. The power law spectrum 
follows from the results of Fermi first order acceleration as- 
sumed for individual accelerators (e.g., Bell 1978, see also 
Ellison et al. 2004). The normalization is obtained from the 
total power transferred by supernovae into the CRs kinetic 
energy within a given volume 



^C inj \-^p, kin j-C'p, kin'2-C'p, kin 



-KEZ p ±\j(-p + 2) 



p, kin, mil 



Y^niPK-i/V, 



(3) 



where, for simplicity, p/ 2 was assumed, and we used the 
fact that i?p, kin, min <S -E Pi kin, max in the second equality. 
The last part of Eq. © is simply the total power using IZi 
(y] - IZi — TV) as the rate of supernova explosions in the star 
forming region of volume V that release a fraction r\i of its 
power V into relativistic cosmic rays. We note that nonlin- 
ear theories of CR acceleration predict deviations from the 
power law injection spectrum, especially at low energies and 
near the high energy cutoff (e.g., Berezhko & Ellison 1999; 
Ellison et al. 2000). In addition, Ptuskin & Zirakashvili 2005 
argue that at any given time, only particles near an evolving 
cutoff energy actually escape the supernova remnant, and 
only the time-averaged injection spectrum is a power law 
(see also Fujita et al. 2009; 2010). Thus, the spectrum and 
time evolution of CR accelerators may be more complicated 
than we consider here. 

The injection spectrum of the primary particles (Q(E), 



from which one derives the distribution N(E) and then e.g., 
the gamma-ray emission) is assumed to be directly normal- 
ized by the total supernovae explosions in the galaxy, and to 
have a power-law form with fixed index, irrespective of po- 
sition within the starburst. The common approach is then 
equivalent to assuming that at any point of the starburst 
region, there is a single accelerator inputing cosmic rays 
at a reduced rate, one that is compatible with the overall 
volume-average given by the total number of accelerators. 
How realistic are all these assumptions? 

There are direct observational hints that the luminos- 
ity in gamma-rays is correlated with the star formation rate 
(from now on, SFR) or supernova frequency. The correla- 
tion may be slightly non-linear, following a power-law with 
a hard index, L 7 ~ SFR 6 , with b ~ 1-1.4 (Abdo et al. 2010, 
2011). Due to the current scarcity of starbursts measured in 
gamma-rays, this is, however, not yet proven. Strictly speak- 
ing, this relation is in fact impossible to predict in the models 
presented in the literature so far, since the SFR goes into the 
normalization of the primary cosmic ray spectrum, which is 
assumed. Then, it is natural to expect that to first order, the 
same correlation is maintained in the hadronic gamma-ray 
luminosity generated by the former population. Plots such 
as L 7 vs. SFR (or similarly, enhancement of cosmic-rays vs 
SFR) have not been derived (to our knowledge) from first 
principles yet. 

We also note that the problem of calorimetric (essen- 
tially, full conversion of the proton kinetic energy into sec- 
ondary particles like gamma-rays) or non-calorimetric be- 
havior also depends on the validity of the assumption of 
the primary spectrum, emphasizing that an understanding 
on how it is built up from individual contributors could be 
useful. 

The aim of this work is then to assess some of the 
formerly-commented common assumptions made in the de- 
scription of the injected cosmic-ray spectrum of starbursts 
regions, by accounting for the individual contributions of ac- 
celerators residing in the starburst volume over the starburst 
age. 



2 NUMERICAL APPROACH AND RESULTS 

As stated, instead of considering a total, volume-averaged, 
cosmic-ray injection, we shall obtain the individual contri- 
bution of each accelerator to the cosmic-ray intensity at a 
given time. The cosmic ray intensity contribution of each 
one of the accelerators will thus be given by 



ME, 



*> = S'' 



(4) 



where f(E, r, t) is the distribution function of protons at 
time t and distance r from the source which satisfies the 
radial-temporal-energy dependent diffusion equation 



(df/dt) = (D(E)/r 2 ){d/dr)r 2 {df/dr) + 

(d/dE) (Pf) + , 



(5) 



(we choose spherical coordinates, with R being the radial 
distance from a given accelerator, and P representing the 
energy losses). We should note the differences between the 
simplified Eq. (0, for the whole starburst, and Eq. ((5j. The 
latter is written for a single accelerator and thus, Q is the 
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Figure 1. Concept for the computation of individual contribu- 
tions of the cosmic-ray density in star forming regions. 



injection of cosmic-rays of that accelerator only and is ob- 
viously more constrained. Note too that Eq. ([5]) refers to an 
E-, R-, and t— dependent problem, whereas after the sim- 
plifications over Eq. {1} that we mentioned above, it is ho- 
mogeneous and in steady state, and thus only E— dependent. 
But we note a cavaet too: in order to work with analytically 
tractable solutions of Eq. J5J, see below, we ignored the ad- 
vective escape term (i.e., N(E)/r(E)); thus, we do not con- 
sider here outflows that may remove cosmic rays from the 
starburst regions. Cosmic rays are transported here by dif- 
fusion, which is appropriate in the calorimetric limit only 
(see, e.g., Loeb & Waxman 2006, Thompson et al. 2007).. 

We assume that the proton energy loss is due to nuclear 
interactions. The nuclear loss rate is P nuc = E/t vv , with 
T pp — (n p cna pp )~ is the timescale for the corresponding 
nuclear loss, k ~ 0.45 is the inelasticity of the interaction, 
and (Jpp ~ 33 nib is the cross section (e.g., Stecker 1971, 
Gaisser 1990); a pp is essentially constant for the energies 
of interest above GeV. Thus, following e.g., Aharonian & 
Atoyan (1996), t pp ~ 7x 10 7 (n/cm" 3 ) -1 yr. The solution to 
the diffusion equation for an arbitrary energy loss term, dif- 
fusion coefficient, and impulsive injection spectrum /mj(-E), 
such that Q(E,r,t) — Nofi n j(E)Sf5(t), for the particular 
case in which D(E) oc E and /i n j oc E~ a , is (see Aharo- 
nian & Atoyan 1996 and references therein) 

f(E,r,t)~(N E- a /iv 3/2 R 3 dit ) 

X exp [-{a - l)t/r pp - (R/R di{ ) 2 ] , 

where 



(6) 



R diS = 2JD(E)t 



^exp(t5/r pp ) - 
(tS/Tpp) 



(7) 



stands for the radius of the sphere up to which the particles 
of energy E have time to propagate after their injection. In 
case of continuous injection of accelerated particles, given 
by Q(E, t) — QoE~ a T(t), the previous solution needs to 
be convolved with the function T(t — t') in the time interval 
^ i! ^ t. If the source is described by a Heavside function, 
T(t) - Q(t) the solution is (Atoyan et al. 1995) 

f(E, r, t) = (Q E~ a /47vD(E)r)(2/V^) 



All previous equations in this section are thus valid for a 
single accelerator, e.g., a single supernova remnant or pul- 
sar wind nebula, injecting either impulsively or continuously 
cosmic rays into the starburst. We now consider that a star- 
burst is a collection of such accelerators, which are assumed 
to appear in the starburst volume at a given rate, e.g., equal 
to the SN rate (which is in turn related to the SFR), and 
which could be constant (as is usually assumed) or not (e.g., 
the starburst phenomenon could have ended, or it could be a 
short burst of star formation rather than a continuous one). 
Thus, one can consider the evolution in time of the injection 
spectrum of the whole starburst by computing each acceler- 
ator's contribution using the solution to Eq. (0), and solving 
it as many times as there are injectors in the volume. 

This concept is depicted in Figure [T] Each plane repre- 
sents a given past time in the evolution of the starburst (the 
current time is £3 in that figure). Each test region of the star- 
burst today receives the cosmic-ray contribution of each of 
the accelerators injecting cosmic-rays in its past, for which 
we individually consider the propagation from the time of 
its injection up to the test region. Summing up all contri- 
butions, we build up the spectrum at that point. Compar- 
ing different test points we can check for homogeneity, and 
study, for instance, how many and which are the accelera- 
tors contributing mostly to the total spectrum. One conse- 
quence of such an approach to compute the total spectrum 
is a direct relation between the accelerator appearance rate 
and the cosmic-ray intensity in a given position. By slicing 
in time and space bins from that test-position, the number 
of accelerators in each bin will, under a homogeneous and 
isotropic distribution, be directly proportional to the accel- 
erator appearance rate. (Each of the panels in Figure [l] will 
have as many accelerators as determined by the appearance 
rate assumed.) Such proportionality would be apparent in 
the cosmic-ray intensity, since, for sufficiently small time and 
space bins, the differences among the contributions of accel- 
erators that were born in that bin will be negligible; this 
would lead to a linear relationship between cosmic-ray in- 
tensity and accelerators appearance rate. We have checked 
that this is indeed the case in our numerical simulations. 



2.1 Cosmic-ray injection from a single accelerator 

As an example of the results, we consider first the injection 
of a single accelerator. We show in Figure [2] the cosmic- 
ray intensity contribution for an injection that occurred at 
different times in the starburst history, and at different dis- 
tances from the test point. We consider here (as we do be- 
low) that the diffusion coefficient at 10 GeV is 10 26 cm 2 s _1 
and 8 = 0.5 in a medium of density n — 35 cm~ 3 ; we have 
also explored other parameter values. A low value of diffu- 
sion coefficient (slow diffusion) is expected in dense molec- 
ular regions, as an starburst (see e.g., the cases of dense 
molecular clouds in our Galaxy in the works by Ormes et al. 
1988, Gabici 2010, or Torres et al. 2010, and other references 
therein). With the selected value of interstellar particle den- 
sity, the total molecular mass in a typical starburst sphere 
of 300 pc is ~ 10 8 M©. Times considered are 10 4 , 10 5 and 
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Figure 2. Cosmic-ray intensity contribution for injection at dif- 
ferent times in the starburst history, and at different distances 
from the test point. Times considered are, from top to bottom, 
10 4 , 10 5 and 5 X 10 5 years. 



5 x 10 years and distances considered are 10, 50 and 100 
pc from the test point. 

Figure [2] can be understood assuming that tS <C 
T pp , and thus, from Eq. and using the adopted 
parameters, one gets the diffusion radius as i?dif ~ 
12 pc (E/10 GeV) 1/4 (t/10 5 yr) 1/2 . We see that only high- 
energy particles diffuse away fast enough to reach distant 
regions. We also note that the ages of the accelerators we 
(purposefully) choose are reflected in low individual contri- 
butions to the cosmic-ray intensity. They have, however, a 
varying E— dependent importance of which injectors are the 
dominant contributors. 

This effect is more clearly seen in Figure 3, in which we 



Figure 3. Ratio of the cosmic-ray intensity as a function of en- 
ergy between a young and an old accelerator; i.e. one that injected 
particles at the earlier times of 10 4 and 10 5 years (top), and 10 4 
and 5 X 10 s years (bottom), for the three selected distances to 
the test point. 



plot the ratio of the cosmic-ray intensity as a function of en- 
ergy, between an accelerator that injected particles 10 4 and 
10 years ago (top), and between 10 and 5 x 10 years ago 
(bottom). At low energies and distances, the contribution of 
the older accelerator dominates that of the younger one. At 
10 GeV, for instance, the diffusion radius of the 10 4 year old 
accelerator is ~ 4 pc; whereas it is ~ 12 pc for the 10 5 year 
old one. At 10 pc then, the older accelerator dominates. At 
higher energies, the diffusion radius of the 10 4 year old ac- 
celerator is larger than the test point considered and starts 
dominating. These results are consistent with the study by 
Aharonian & Atoyan (1996), and are just put here in the 
context of injection ages that are much larger than those 
typically considered for a single source (e.g., when study- 
ing the interaction of such cosmic-rays with the surrounding 
molecular gas, like in mid-age SNRs). 



2.2 Homogeneous starbursts at constant 
explosion rate 

Taking into account the individual contributions, we start by 
supposing that the starburst is a sphere (e.g., with a radius 
of 300 pc) and has experienced a constant supernova explo- 
sion rate (e.g., of 0.1 SN year^ 1 ) during the recent history 
(e.g., during the last 1 x 10 6 years). With the previous values, 
we have a volume- average number of explosions (generically 
called accelerators or injectors) equal to ~ 9 x 10 -4 pc -3 . In 
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accord with the approach described above, the total cosmic 
ray distribution is 



AT 
2 , -7-fn(E,d n ,t n ) 



(9) 



where f n is the individual cosmic-ray contribution. In this 
example, then, we numerically consider 10 J individual solu- 
tions to Eq. (J5]). For this purpose we developed a numeri- 
cal scheme that solves the diffusion-loss equation for each 
injector, under the assumption that the cosmic-ray source 
is impulsive if the time past since its appearance is larger 
than the adopted regime break of 10 4 years and continuous 
otherwise, and sums them all up. We randomize via Monte 
Carlo the position and time of appearance of each accel- 
erator (with the specified mean explosion rate), as well as 
the individual parameters, in particular, the energy input is 
randomized between (0.1 — 1) x 10 50 erg for an impulsive 
case, and between (0.1 — 1) x 10 37 erg s _1 for a continuous 
case. We again assume that the diffusion coefficient at 10 
GeV is 10 26 cm 2 s _1 and S = 0.5 in a medium of density 
n = 35 cm -3 . Each source spectrum was taken to have a 
slope of 2.2, but we also considered the case where the slope 
was randomly selected in an interval around this value. 

Our results verify the approximate isotropy of the 
cosmic-ray distribution. Unless we are interested in an analy- 
sis test point that is too close to a recent accelerator location, 
we find that the cosmic-ray distribution is homogeneous. For 
example, in the top panel of Figure [4] we show an example 
of the results for antipodal positions in azimuth and zenith 
angle, at a distance of 100 pc from the center of the star- 
burst. Our results also show that the cosmic-ray intensity 
is a function of starburst-center distance, being smaller the 
farther away from it we locate the test region. This is the 
result of a finite-size starburst sphere, and appears as a re- 
flection of the fact that test regions closer to the limit of 
the star-formation environment lack the contributions of in- 
dividual accelerators located beyond this distance. This is 
shown in the bottom panel of Figure [3] for test positions lo- 
cated at different distances from the center of the starburst. 
In configurations such as this, namely a starburst sphere 
with a limited size, the difference between the cosmic-ray 
spectrum at the center and its boundary can reach up to a 
factor 2. This is usually neglected in the models published 
in the literature (see the introduction). 

We compute the cosmic-ray intensity also as a function 
of time from the beginning of the burst, and an example 
is shown in Figure [5] There we show the cosmic-ray inten- 
sity in the starburst region at 50 pc from the center at dif- 
ferent times, from 5000 (bottom curve) to 10 7 years. We 
note that as long as the starburst is ongoing, the cosmic- 
ray intensity increases at all energies. Although the rate of 
increase of the cosmic-ray intensity slows down after a few 
hundred thousand years, it only attains steady-state after 
a few million years. For ~ 10 7 years, the cosmic ray inten- 
sity in Figure [5] would overlap with the curve corresponding 
to 10 6 years. It is then important to know the age of the 
starburst when detailed models are constructed, at least for 
ages up to T pp ~ 10 6 years. The latter actually depends on 
how dense the starburst is, in other words, on how the age 
of the starburst compares with t pp . 
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Figure 4. Top: Cosmic-ray intensity for antipodal positions in 
azimuth and zenith angle, all at a distance of 100 pc from the cen- 
ter of the starburst. Errors result from a small number of simula- 
tions (6 realizations only), and are shown merely for illustration. 
The larger the number of simulations, the closer all curves are 
to one another. For the sake of clarity of presentation, the (ordi- 
nate) flux was multiplied by the factor E 22 . Bottom: Cosmic-ray 
intensity for test region positions located at different distances 
from the center of the starburst. 
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Figure 5. Cosmic-ray intensity in the starburst region, at 50 pc 
from the center, and at different times, from 5000 (bottom curve) 
to 10 7 years. The curve for 10 7 years overlaps the one for 10 6 
years. Errors are derived as the standard deviation for each case 
after hundreds of MC realizations. 
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Figure 6. Top: Example of cosmic-ray intensity contributions as 
a function of distance (around 50 pc from the test position, lo- 
cated at 100 pc from the center of the starburst). The blue points 
show accelerators that contribute more than 10% of the total flux, 
in green, those contributing between 1 and 10%, and in red those 
contributing between 0.1 and 1%. Bottom: Contributions of in- 
dividual accelerators surrounding a test-position at 100 pc from 
the center; black dots represent contributions of less than 0.1 %. 



In our current setting and for our assumed explosion 
rate, we have about ~ 500 accelerators exploding over the 
starburst age at positions within 50 pc of any location of in- 
terest (except of course those positions located at the bound- 
ary of the starburst). These accelerators provide more than 
97% of the cosmic ray intensity at each point. This can be 
seen in the example of Figure [6] where we plot the percent- 
age of the total (integral) cosmic-ray intensity contributed 
by individual accelerators as a function of their distance to 
the test point. The contributions of individual accelerators 
surrounding a test position at 100 pc from the center are 
shown in the bottom panel of Figure \fi\ 



2.3 Varying the diffusion coefficient 

Since the diffusion coefficient is actually unknown, we have 
also considered the effect of a faster diffusion (a larger 
DioGcV , leading to a shorter timescale for diffusion of parti- 
cles). As expected, the larger D is, the smaller is the cosmic- 
ray population at higher energies. This is exemplified at a 
fixed distance (50 pc) from the center of the starburst in 
Figure [7] We also plot in the bottom panel of that Figure 
the cosmic-ray intensity for test positions located at differ- 
ent distances from the center of the starburst, as in Figure|4] 
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Figure 7. Top: Example of cosmic-ray intensity at 50 pc from 
the center as a function of D. Bottom: Cosmic-ray intensity for 
test positions located at different distances from the center of the 
starburst, as in Figure|4]but for -DioGcV = 10 27 cm 2 s _1 . 



but for DioGcV = 10 27 cm 2 s _1 . These figures lead to the fol- 
lowing conclusion: It is interesting to note that the discovery 
of a hard spectrum up to the highest energies would favor a 
small diffusion coefficient in starburst galaxies. If one consid- 
ers that gamma-rays are ~10 times less energetic than the 
protons that generated them, differences in the gamma-ray 
spectrum at say, 1-10 TeV should be clear despite possi- 
ble degeneracies introduced by differences in modeling the 
injection or the size of the emitting region. Generally, the 
flatter the spectrum is at the highest energies, the better 
is the case for a small D. However, it is difficult to deter- 
mine such changes from analysis of current measurements. 
Figure [7| shows that for Din GeV = 10 27 cm 2 s _1 , the E 22 J P 
flux drops by a factor of 4 over an energy range of 1000, 
between E p — 10 GeV and 10 TeV (roughly corresponding 
to E~f ~ 1 GeV and 1 TeV). That translates to a spectral 
steepening of 0.2. This measurement is particularly suitable 
for the forthcoming Cherenkov Telescope Array (Actis et al. 
2011), which should provide a detailed gamma-ray spectrum 
of nearby starburst galaxies from a few tens of GeV up to 
100 TeV. 



2.4 Stability of results 

We have tested many realizations of the set of positions 
where the injectors are located along the lifetime of the star- 
burst, and the only visible difference is set by the proximity 
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Figure 8. Cosmic-ray intensity at 50 pc from the starburst center 
for an accelerator appearance rate of 0.01 (top) and 0.1 (bottom) 
and derived standard deviation after thousands of MC realiza- 
tions. 



to the test point of the nearest injector. This produces differ- 
ences especially at low energies, given the results commented 
in the previous section, but not an appreciable difference at 
energies beyond 100 GeV. At high energies, since the dif- 
fusion scales grows, the contribution is about the same ir- 
respective of the distribution (high-energy protons diffuse 
significantly away from accelerators). As an example of the 
stability of the results obtained out of a given realization 
we plot the cosmic-ray intensity at 50 pc from the starburst 
center for an accelerator appearance rate of 0.01 and 0.1 (see 
Figure [5}. The errors are the derived standard deviation af- 
ter thousands of MC realizations. Points at higher energies 
show barely any error. 

We ran two additional tests. The first one is on the 
influence of a fixed value of the slope of the cosmic-ray in- 
jection for all accelerators instead of a variable one, which 
we have chosen in the range 2.0 - 2.4 (see, e.g. Bell et al. 
2011). That is, in the second case each injector has a dif- 
ferent (randomly assigned) slope on its cosmic-ray spectral 
injection. The results are depicted in Figure [5] where we see 
that the change is visible especially at high energies, where 
it can reach up to ~ 40 % due to the influence of harder 
injectors. 

We have also tested the influence of a different cut be- 
tween the impulsive and continuous assumptions for the ac- 
celerators (i.e., a different regime break), by studying the 
stability of the results for breaks at 50000 and 100000 years. 
The second panel of Figure [9] shows the ratio for the cosmic- 



Figure 9. Ratio of cosmic-ray intensity for configurations al- 
lowing a variable [2.0-2.4] and fixed [2.2] slope of the injection 
spectra, and its derived standard deviation after thousands of 
MC realizations of the accelerators positions and times of their 
appearance. Botton: Ratio for cosmic-ray intensity with different 
regime breaks. In both cases the test point was chosen at 50 pc 
from the center of the starburst, as an example. 



ray intensity with different regime breaks with that consid- 
ered earlier at a test point 50 pc away from the starburst 
center. The larger the regime break, the more accelerators 
are considered continuous injectors. This explains why the 
cosmic-ray intensity is smaller the larger the break point: 
the higher the energy, the shorter is the time a particle has 
to reach the test point in the case of a continuous accelera- 
tor, producing a difference of up to a few percent in the case 
of a 50000 years break. 



2.5 Terminated star-formation burst 

To exemplify the cosmic-ray decay, we analyze the situation 
in which the star formation is exhausted in the whole star- 
burst. Figure [10] shows how fast the cosmic-ray spectrum 
dies out when the star formation stops as a function of time 
from the end of the burst. It is interesting to note that the 
cosmic-ray spectrum, even in cases in which the last accel- 
erator appears 5 x 10 5 years ago, is only a factor of 2 less 
than what it would be under a continuous star forming pro- 
cess. The timescale for decay in this model is of order ~ t pp , 
as might be expected if pion production is the only loss. 
The decay of the cosmic ray spectrum is slow as shown in 
the bottom panel of Figure 1101 We obtained similar results 
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Figure 10. Top: Cosmic-ray spectrum today for different times 
at which the star formation is exhausted. Bottom: decay of the 
cosmic-ray intensity as a function of time and energy. 



also for the case in which there is a remaining, residual star 
formation only in the central part of the original region. 
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Figure 11. Model for cosmic-ray distribution in NGC 253. Shown 
in the upper panel are proton spectra at high energies, and at 
different distances from the galaxy center, along the plane. The 
bottom panel shows the ratio between the cosmic-ray distribution 
at 50 pc in the ISB and the one at different distances, within the 
central core and along the disk, as a function of energy. Two 
different accelerator appearance rates are considered. 



3 THE CASE OF NGC 253 

We now consider a core-enhanced star-formation rate distri- 
bution where the starburst geometry is divided into a central 
sphere (or innermost starburst), in which the star forma- 
tion is maximal, and a surrounding disk. Such a situation 
is typical in starbursts galaxies, and has been considered, 
for instance, for NGC 253, where infrared, millimeter, and 
centimeter observations show that the central region dom- 
inates ongoing star formation (e.g., Ulvestad & Antonucci 
1997; Ulvestad 2000). Current estimates of the gas mass in 
the central 20-50 arcsec region range from 2.5 x 10 7 Mq (Har- 
rison et al. 1999) to 4.8 x 10 8 Af Q (Houghton et al. 1997); 
see Bradford et al. (2003), Sorai et al. (2000), and Engel- 
bracht et al. (1998). We also note that a detailed numeri- 
cal treatment of the steady-state distribution of cosmic ray 
electrons and protons in NGC 253 was recently carried out 
by Paglione et al. (1996), Domingo-Santamaria & Torres 
(2005), and Rephaeli, Arieli, & Persic (2010). In the latter 
case, their modified GALPROP code follows the evolution 
of particle energy and spatial distribution functions in the 
context of a diffusion-convection model. No time-dependent 
analysis of the contribution of individual accelerators is in- 
cluded. 

Here, following Domingo-Santamaria & Torres (2005) 



we will assume, in agreement with the mentioned measure- 
ments, that within the central 200 pc (100 pc radius), a disk 
of 70 pc hosts 3 x 1O 7 M0 of molecular mass, uniformly dis- 
tributed, with a density of ~ 600 cm" . In this region, the 
supernova explosion rate is taken as 0.08 year -1 . We have 
also considered a lower value for the accelerator appearance 
in the core (0.02 yr -1 ). This is motivated by the relationship 
between supernova rate and IR luminosity (Tacki et al. 2011) 
and also by the fact that Melo et al. (2002) indicate that 50% 
of the IR luminosity (and 70% of the star-formation) is in 
the starburst and 30% of the IR and star-formation is in the 
surrounding disk. Williams & Bower (2010) also find that 
about half of the GHz radio emission is extended (from the 
disk) and half comes from the starburst core, which, to the 
extent the FIR-radio correlation holds, suggests a roughly 
equal star-formation rates in the starburst and the surround- 
ing disk. 

Additional mass with an average density of ~ 50 cm -3 
is assumed to be in the central kpc outside the innermost 
region, but subject to a smaller injector rate of 0.01 year -1 , 
or ~ 10% of that found in the most powerful nucleus (Ulves- 
tad 2000). We assume that each injector accelerates cosmic- 
rays with a spectral slope between 2.2 and 2.4, randomly 
assigned. 

Figure [Til shows the result for the cosmic-ray distribu- 
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tion at high energies, and at different distances from the 
center of NGC 253, within and outside the innermost star- 
burst region with the highest accelerator appearance rate. 
Results for different accelerators appearance rate in the in- 
nermost starburst are shown. Each value of N p is obtained 
by summing up all contributions coming from accelerators 
located in the innermost starburst as well as those located 
in the disc. The top panel also shows the average cosmic- 
ray distribution obtained by Domingo-Santamaria & Torres 
(2005), based on the treatment described in the introduc- 
tion. The impact of a core-enhanced rate distribution is ev- 
ident, with a significant change, of more than an order of 
magnitude at 1 TeV, in the cosmic-ray density when the 
distance is increased from 50 pc to 200 pc from the center 
(i.e., from within to outside the innermost nucleus). There 
is in addition a significant change in the cosmic-ray distri- 
bution further outside of the nucleus, the plot shows it for 
distances from 200 to 400 pc along the surrounding disc. 
The bottom panel of Figure [TT1 shows the ratio between the 
cosmic-ray distribution in the sphere and the one at differ- 
ent distances from the center as a function of energy. It can 
seen be again that large values of the ratio are obtained; 
even for a smaller difference between the rates of the (in- 
nermost starburst) ISB and that of the surrounding disc. 
Thus, one-zone models are a rough approximations to this 
situation. 



4 CONCLUDING REMARKS 

The results presented in this paper are based on a step- 
by-step approach for the computation of the build-up of 
the cosmic-ray spectrum in star forming environments: The 
respective contribution of each accelerator injecting cosmic 
rays over the history of the star forming region is individ- 
ually accounted for. The total spectrum is then obtained 
by summing up all individual contributions, taking into ac- 
count an energy, time, and space-dependent propagation of 
the cosmic-rays in the starburst environment. Even with the 
caveat of not yet considering outflows, the total spectrum 
shows two features that are found from gamma-ray observa- 
tions of starburst galaxies, 1) that the average photon index 
of the current cosmic-ray spectrum in starbursts is hard, 
even though the cosmic ray spectrum of the individual in- 
jectors suffer significant steepening due to their propagation, 
and 2) put in context the fact that the integrated cosmic-ray 
spectrum (and its gamma-ray yield) scale linearly with the 
star formation rate. The latter can also be considered an ob- 
served fact (see Abdo et al. 2011) given that the gamma-ray 
luminosity of nearby star forming galaxies seem to scale lin- 
early with the radio continuum and infrared luminosity, both 
established tracers of star formation. Our simulations shows 
that the discovery of a hard spectrum up to the highest en- 
ergies would favor a small diffusion coefficient. Our results 
also prove that in the calorimetric limit (with no cosmic-ray 
advection), homogeneity is reached (typically within 20%) 
across the whole starburst region. However, values of center- 
to-edge intensity ratios can amount to a factor of several. 
Differences between the common homogeneous assumption 
for the cosmic-ray distribution and our models are larger in 
the case of two-zone geometries, such as a central nucleus 
with a surrounding disc. We have also found that the de- 



cay of the cosmic-ray density after the starburst process is 
terminated is slow, and even approximately 1 Myr after the 
burst ends (for a gas density of 35 cm~ 3 ) it may still be 
within an order of magnitude of its peak value. 
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